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Abstract 



Interacting surfaces are frequently found in mechanical systems and components. A lubri- 
cant is often added between the surfaces to separate them from mechanical contact in order 
to increase life and performance of the contacting surfaces. In this work various aspects of 
hydrodynamic lubrication are investigated theoretically. This is where interacting surfaces 
are completely separated by a fluid film which is often the desired operating condition of 
machine components when wear and friction is to be reduced. Different flow regimes can 
be identified within the scope of hydrodynamic lubrication. 

If the surfaces are separated by a thick fluid film the influence from surface asperities is 
small and the surfaces can be treated as smooth. If the rate of change in film thickness with 
respect to the spatial directions is significantly large and if the flow velocity or Reynolds 
number is large, the ordinary fluid mechanical approach treating viscous flow with Com- 
putational Fluid Dynamics (CFD) has to be used. CFD is used to investigate influence 
from the use of an artificial microscopic surface pattern on one of the two interacting sur- 
faces. The influence from the pattern is isolated from any other pressure generating effects 
by keeping the interacting surfaces parallel. Results are shown for different shapes of the 
micro-pattern. If the Reynolds number decreases, the system enters a regime called Stokes 
flow where the inertia effects are neglected. The full CFD approach is compared with the 
Stokes for various physical and geometrical cases. 

If the change in film thickness is small in the spatial directions, the thin film approxi- 
mation is applicable and the full momentum equations describing fluid flow together with 
the mass continuity equation can be reduced to the Reynolds equation. Depending on 
boundary conditions, low pressures can occur at location of expanding fluid gap leading to 
tensile stress applied to the lubricant. However, a real liquid lubricant can only resist small 
tensile stresses until it cavitates into a mixture of gas and liquid. This often happens close 
to atmospheric pressure due to contamination and dissolved air into the liquid and oc- 
curs at higher pressures than the actual vaporization. To avoid pressures reaching too low 
levels, a general cavitation algorithm applied to the Reynolds equation is presented that 
accommodates for an arbitrary density-pressure relation. It is now possible to model the 
compressibility of the lubricant in such a way that the density-pressure relation is realistic 
through out the contact. The algorithm preserves mass continuity which is of importance 
when inter-asperity cavitation of rough surfaces is considered. 

For small film thicknesses the surface roughness becomes important in the perfor- 
mance of the lubricated contact. Even the smoothest of real surfaces is rough at a micro- 
scopic level and will influence the contact condition. The Reynolds equation still applies 
since the heights of the surface asperities are small compared to the spatial elongation. 
Treatment of the roughness of a real surface in a deterministic fashion is however beyond 
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the scope of today’s computers. Therefore other approaches need to be employed in or- 
der to take the surface roughness into account. In this work a homogenization method is 
used where the governing equation of the flow condition is formulated with a two-scale 
expansion, the global geometry and the roughness. Solutions are achieved for the limit 
of the roughness wavelength, e - - 0 and the method renders a possibility to treat the two 
scales separately. A method to generate dimensionless flow factors compensating for the 
surface roughness is developed. The flow factors, once solved for a particular surface, can 
be used to compensate for the surface roughness in any smooth global problem for any 
film thickness. 
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Nomenclature 



0* Couette flow factor 

Ot V max! h 

0s Poiseuille flow factor 

h{x) Nominal fluid film clearance 

p Dimensionless pressure, p/p c 

x Dimensionless x-coordinate, (6r|Mx)/(p(z mflA: ) 2 ) 

P Bulk modulus 

X Solution to the local problem 

Ax Local grid size 

r| Dynamic viscosity 

A 6t|m/P 

L x , L y Typical length scales in x-, v-direction 

u(x,y,z) Velocity field 

\|/ Solution to the local problem 

p Mass density 

p c Density at atmospheric pressure 

a Standard deviation 

9 Dimensionless density, p/p c 

r Spatial coordinate or surface topography height 

r max Maximum surface asperity height 

£ Wavelength measure of a local surface roughness scale 

\ x/e 



m 

Pa 

m 

Ns/m 2 

m 

m/s 

kg/m 3 

kg/m 3 

m 

m 



IX 



Awall 


Wall area, L x x L z 


2 

m 


d 


Groove depth 


m 


d + 


d/Ly 




F x 


Tangential wall force 


N 


F+ 

1 X 


Dimensionless tangential wall force, 100 * F x L? y / (, A wa ur\uL x ) 




Fy,W 


Load carrying capacity 


N 


F y + 


Dimensionless normal wall force, 100 * F y Ly / (A wa nr\uL x ) 




G 


Dimensionless surface topography height, —z/zmax 




8 


Switch function, 0 if 0 < 1 and 1 if 0 > 1 




H 


Dimensionless film thickness, h/h 




h(x) 


Film thickness, 


m 


h 0 


Minimal fluid film clearance 


m 


K 


Kurtosis of surface topography 




8 V- A- /- 


Fluid domain length in x-, y- and z-direction 


m 


N 


Number of grid nodes for the cell problem 




P 


Pressure 


Pa 


P + 


Dimensionless pressure, 100 * pL y / (r\uL x ) 




Pc 


Atmospheric pressure 


Pa 


Pe 


Peclet number, pwAx/r| 




Ra 


Average surface roughness 


/j m 


R, 


RMS surface roughness 


jjm 


Re 


Reynolds number, puL y /r[ 




SK 


Skewness of surface topography 




u, v, xv 


Velocity in x-, y— and z-direction 


m/s 


ui(x,y,z) 


Velocity field tensor 


m/s 


XV 


Groove width 


m 


XV+ 


xv / L x 




X, y,z 


Spatial coordinate directions 


m 


x + 


Dimensionless x-coordinate, 2 x/L x — 1 




xd 


Groove displacement 


m 


xd + 


xd/xv 





Chapter 1 

Introduction 



The widespread range of mechanical applications with interacting surfaces and an in- 
creased industrial competition for improved performance and more cost effective products 
have lead to an increased attention to research in the science of friction, wear and lubri- 
cation, i.e. tribology. The rapid improvements of computers allow for more advanced 
tribological problems to be solved in a numerical fashion. This has led to an increased 
potential of numerical research and virtual experiments as a supplement to experimental 
approaches. Theoretical research is beneficial in terms of reducing expensive tribological 
testing and opens opportunities for more detailed investigations of the influence of spe- 
cific physical parameters in order to increase the understanding of the physical problem 
considered. 

In this thesis theoretical research is carried out to cover various aspects of sliding 
lubrication. A vast number of mechanical applications contain interfaces with sliding in- 
teracting surfaces. These interfaces operate under various loading conditions and relative 
sliding speeds. To reduce friction and minimize wear, the interacting surfaces are often 
separated by a thin layer of lubricant. Depending on the surface geometries, relative speed, 
surface topography and applied load among other things, the surfaces could be partly or 
completely separated from each other. Surfaces must form a converging gap or wedge to 
generate an effect of separation. When the lubricating fluid is dragged into the converging 
gap because of the relative motion of the surfaces, pressure will rise due to fluid compres- 
sion because of the decrease in film thickness, producing a hydrodynamic load carrying 
capacity of the lubricated contact. 

There are several lubrication regimes describing different operating conditions of the 
contact. An instructive illustration is shown in Fig. 1.1. The surface mean separation h 
scaled with the standard deviation of the surface roughness a is plotted as a function of 
the lubricant viscosity and relative velocity divided by the applied load on the surfaces. 
For highly loaded lubricated surfaces at a small relative sliding speed, the regime is called 
boundary lubrication. In this regime all of the load is carried by the asperities of the sur- 
faces. If the speed increases or load decreases the film thickness increases and the regime 
becomes the mixed lubrication. At this regime the load is carried by both the lubricat- 
ing fluid and asperity contact. Next to the mixed lubrication regime is the hydrodynamic 
lubrication (HL) which is considered in this thesis. Within this regime the surfaces are 
completely separated by the fluid film and load is carried by hydrodynamic action. 

Contacts operating in the hydrodynamic regime are found in various machine elements. 
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CHAPTER 1. INTRODUCTION 




Figure 1.1: Different lubrication regimes. 



Typical examples are the conformal hydrodynamic slider bearings, for example journal 
and thrust bearings. They are often designed to carry high loads and are typically found 
in engine crank shafts and turbines. Fig. 1.2 shows a thrust bearing having six pads with 




Figure 1 .2: A typical tilting pad thrust bearing. 

self-adjusting tilt angle. The pads of the bearing are loaded against a collar on a rotating 
shaft. When the bearing is operating the pads are self-adjusted to form a converging gap. 
The lubricant is pulled into the interface to produce pressure and a load carrying capacity. 

Bearing design is complex and involves optimizing several parameters such as the 
bearing geometry, lubricant properties and flow rate. The surface topography affects the 
performance of the bearing and is thus also an important design parameter. Simplistic de- 
sign tools are available in the engineering community but their solutions are often limited 
in accuracy. More general comprehensive equations taking into account surface topogra- 
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phy are often difficult and time consuming to solve which limits their use in industry. For 
this reason theoretical research is important to gain more understanding of the effects of 
surface topography on a lubricated contact performance and to find more effective tools of 
bearing design. 



Chapter 2 

Hydrodynamic Lubrication 



To theoretically describe the nature of a hydrodynamically lubricated contact the physics 
needs to be mathematically interpreted. The complexity of nature implies the need of a 
simplified interpretations. This is a process of deciding which physical parameters could 
be neglected and which must be retained in order to solve the problem in a way that is 
satisfactory. The simplifications of the problem also open an opportunity to isolate certain 
physical parameters and to study their specific effects on the system. 

In the subject of hydrodynamic lubrication some important output parameters are fric- 
tion force, pressure and load carrying capacity of the surfaces, film thickness, temperature 
and the flow field. Friction force and load carrying capacity are properties that can directly 
be related to the performance of the bearing. Most often a small coefficient of friction is 
desired in a bearing. 

In this thesis, is the main focus is pressure and load carrying capacity. But friction force 
and velocity field are also considered. To treat these properties various fluid mechanical 
problems are considered by solving equations governing different regimes of fluid flow. 
The equations considered in this work are derived and discussed in Appendix A. 

The physical scales in a hydrodynamic contact are of importance. Consider Fig. 2.1 
which is a schematic view of a lubricated contact of two wavy surfaces. The lower surface 
is sliding with a velocity u dragging liquid lubricant (White) with viscosity q and density p 
into the contact. L x and L y are typical length scales for the global geometry irregularities. 
When modeling such a contact there are several important factors to consider. The values 
of T|, p, u, h and the ratio L y l L x play important roles. Also the surface topography visual- 
ized in the enlarged elliptic area is an important factor to consider. It should be noted that 
the scales of the surface topography are highly exaggerated. The local topography ratio 
of scales is typically of the order of L y / L x ~ 10 2 . All these properties determine which 
governing equations can be used to simulate the fluid mechanics. Let us categorize three 
types of equations valid in different sub regimes in hydrodynamic lubrication, namely: 

1. Navier-Stokes equations 

2. Stokes equations 

3. Reynolds equation 

All these sub regimes are affected by the surface topography if the film thickness h is small. 
These equations are derived and discussed in Appendix A. The Navier-Stokes equations 
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CHAPTER 2. HYDRODYNAMIC LUBRICATION 




Figure 2.1 : Different scales of a lubricated contact. 



are the most general and are applicable to all hydrodynamic sub regimes. The Stokes 
equations which are simplifications of the Navier-Stokes are applicable to viscous flow 
with small Reynolds numbers. Reynolds equation is the most simplified and is considered 
in thin film flow with small Reynolds numbers. The Stokes equations are valid whenever 
the Reynolds equation is valid. An attempt to visualize which sub regimes are govern by 
these equations considering pressure as the output variable is shown in Fig. 2.2. It is a map 
with the mean film thickness scaled with the standard deviation of the surface asperities 
on the vertical axis and the Reynolds number Re times the length scale ratio of the surface 
irregularities on the horizontal axis. The gradient color visualizes that the influence of 
surface asperities increases as the ratio h/a decreases. The three different flow equation 
regimes are divided by the dashed lines. It is important to keep in mind that the locations 
of the divider lines are not exact and only reflects the trends to some extent. The work in 
the appended papers A, B and C are categorized in terms of the hydrodynamic regime in 
the map. Fig. 2.2. 

All sub regimes within the map signifies high viscosity flow since the film thickness 
h is small compared to other spatial coordinates. The Navier-Stokes flow regime, see 
Fig. 2.2, signifies large Reynolds numbers and also a rapid change in film thickness with 
respect to the spatial coordinates. Here inertia forces of the fluid are significant. If the 
Reynolds number and the geometrical irregularities L y j L x are decreased, inertia forces 
become negligible and the regime is switched to the Stokes flow sub regime. By decreasing 
the film thickness the sub regime becomes the thin film flow govern by the Reynolds 
equation. If the viscous effects can negligible the Euler equations are valid, however, they 
are not applicable in the field of hydrodynamic lubrication. 
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Figure 2.2: A qualitative map of hydrodynamic fluid sub regimes and the governing equa- 
tions. The locations of the appended papers A, B and C in this map are also shown. 
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CHAPTER 2. HYDRODYNAMIC LUBRICATION 



Chapter 3 

Research Objectives 



It is known that the surface roughness may affect lubricated contact performance in various 
ways. Even if the surfaces are smooth on a macroscopic level, asperities on a microscopic 
level may affect the overall performance in terms of pressure build-up, film thickness and 
friction. It is of importance to optimize surface topography to improve machine perfor- 
mance. This requires better fundamental understanding of contact conditions. 

Deterministic computations of rough tribological contacts are difficult even with mod- 
ern computers due to a large number of degrees of freedom imposed by the roughness. 
Thus, research aiming at finding simplified models that take into account surface rough- 
ness, thermal influence and asperity contact becomes important in order to develop effec- 
tive computational tools for industrial applications. 

The objectives of the present work are as follows: 

1 . To investigate how surface texture can improve sliding contact performance 

2. To investigate the effect of surface topography on thin film contact performance 

3. To develop engineering tools for design of efficient sliding contacts of rough sur- 
faces 



9 



CHAPTER 3. RESEARCH OBJECTIVES 



Chapter 4 

Summary of Appended Papers 



Paper A 

In Paper A the effect of introducing an artificial surface texture or a micro-pattern to one 
of the interacting surfaces is investigated. The pattern is modeled as a pocket and the 
ratio of scales, L y f L x , of the pocket is in the order of 0.1 to 1 which is a large value 
when speaking of thin fluid film flows. Thus, the Navier-Stokes equations have to be used 
since the effects of fluid inertia are large. The influence from the pattern is isolated from 
any other pressure generating effects by keeping the interacting surfaces parallel. This 
approach provides a way to study of how surface texture affects the flow field and pressure 
for different Reynolds numbers and geometries of the pockets that constitutes the surface 
texture. 

Paper B 

It was shown in Paper A that pressure can be low in the diverging part of the pocket if 
the boundary pressure is low. A liquid lubricant cannot withstand tensile strain and will 
cavitate into a mixture of liquid and gas. This happens due to contamination particles 
and gas desolved in the liquid. For this reason the cavitation has to be taken into account 
when treating certain geometries and textures, typically if there is a considerable divergent 
part of the domain. It may also become important to consider the lubricant liquid as 
compressible. In Paper B a cavitation algorithm is developed where an arbitrary lubricant 
compressibility relation in the form p = p (p) can be treated. This algorithm is developed 
for the thin film regime. In this work the influence from the surface asperities is neglected. 
This assumption is valid for values of h/o that are sufficiently large. Fig. 2.2. However, 
with smaller surface separation the influence of surface asperities becomes important. 

Paper C 

In this paper an approach to deal with the influence of surface asperities is considered. 
Fig. 2.2. The thin film regime is considered and it should be noted that the local ratio 
of scales of the surface asperities is typically in the order of L y f L x ~ 1 0 2 which allows 
the approach. Even if the surface is smooth on a macroscopic level the surface roughness 
on a microscopic level will often influence the performance of a bearing in terms of load 



11 



12 



CHAPTER 4. S UMMARY OF APPENDED PAPERS 



carrying capacity. Discreetly resolving the asperities and performing deterministic com- 
putations of the hydrodynamics demands computing powers beyond the scope of today’s 
computers due to the numbers of degrees of freedoms involved. Here a homogenization 
technique is used where the global geometry is separated from the local roughness scale by 
a two-scale asymptotic expansion. The local scale is an arbitrary periodic function with a 
wavelength e — > 0. With this approach a method of computing flow factors, compensating 
for the surface roughness, is developed. The flow factors can be computed for a measured 
surface topography with high resolution and once computed, they can be used as scale 
factors in a coarsely resolved global bearing problem. 



Chapter 5 

Conclusions 



In this thesis, simulation of hydrodynamic lubrication has been considered. Different ap- 
proaches in solving the fluid flow in the contact has been made. 

1 . The fluid mechanics of hydrodynamic lubrication between parallel surfaces is stud- 
ied for isothermal, incompressible and steady 2D conditions. Normally no pressure 
rise in the film would occur when the surfaces are parallel. However, it is shown 
that the introduction of a micro groove on one of the surfaces affects the flow and 
pressure pattern. This gives a net pressure build-up and a load carrying capacity of 
the film. Negligible hydrodynamic effects are achieved when the advective terms 
are excluded from the Navier-Stokes equations, i.e when the Stokes equations are 
considered which proves the significance of inertia for this type of flow condition. 

It is shown that load carrying capacity increases with Reynolds number and the 
width of the groove. A general optimum is achieved when the depth of the groove is 
between 50 to 75 % of the filmthickness for all geometries and groove widths. The 
load carrying capacity can, to some extent, be related to the amount of circulation in 
the groove. A vortex appears at a certain value of groove depth. Near this value the 
maximum load carrying capacity is achieved. 

2. When some part of the fluid gap is expanding pressure can become low and even 
cavitate into a mixture of liquid and gas. For some conditions it also becomes im- 
portant to take the lubricant compressibility into account. To treat this, a generalized 
Reynolds equation accommodating for an arbitrary compressibility relation and cav- 
itation was been derived. It has been shown that different compressibility models 
can produce significantly different results in a hydrodynamic pressure range. Load 
carrying capacity is especially sensitive to the type of compressibility relation used 
if the inlet is starved. Surface roughness close to the inlet has been shown to have 
important effects on load carrying capacity. 

3. For small film thickness the surface roughness becomes important in most cases. 
Due to the amount of degrees of freedom required to resolve the surface roughness 
a homogenization technique for the compressible Reynolds equation was consid- 
ered. With this technique a novel method of computing dimensionless flow factors 
compensating for the effects from arbitrary surface roughness has been developed. 
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By using the flow factors, the hydrodynamics can be solved on a coarsely resolved 
bearing geometry. 



Chapter 6 

Future Work 



Computations of hydrodynamic problems involving surface roughness are complex due to 
the number of degrees of freedom involved in resolving the roughness. The future work is 
to continue studying various effects that surface roughness introduces in lubrication prob- 
lems. The work will be extended to mixed lubrication regime. More physical quantities, 
such as thermal influence will be introduced in the theoretical analysis. Suggested future 
work will also cover developing engineering tools that can be used to predict different 
contact conditions and improve machine performance. These tools should be fast and able 
to cover most lubrication regimes. 
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Appendix A 

Fluid Flow Equations 



The fluid flow in hydrodynamic lubrication is govern by various equations depending on 
the physics to model. These are often described by non-linear partial differential equation 
and solving these equations describing the fluid motion is often complex. Therefore it 
is always desirable to simplify the modeling of fluid flow whenever possible. This is a 
process of determining scales of importance and dropping vanishing terms in the model 
which are not affecting the overall solution. 

The basic governing equations for fluid motion of a viscous fluid are known as the 
Navier-Stokes equations of motion. They are derived by considering the dynamic equilib- 
rium of forces acting at any given region of the fluid and involves conservation of mass 
and momentum. 



A.l Forces 



The forces considered are surface forces, body forces and inertia forces. Body forces 
are proportional to the mass of the fluid and are distributed throughout the mass, e.g. 
gravity force field. Here we will neglect the influence of body force. Surface forces are 
proportional to the extent of the area of the fluid and are exerted on an area element by the 
surroundings through direct contact. The surface forces can be expressed in terms of the 
stress tensor 



T 



ij ~ 



dA 



(A.l) 



where F is the force tensor, A a control area and the first index denotes the components 
normal to the surface and the second components tangential to the surface. In three spatial 
dimensions the stress tensor has nine elements and completely specifies the stress at a 
point anywhere. The stress tensor is symmetric, i.e.. 



Xij — Xjf (A. 2) 

and thus has six independent components for three spatially dimensions. Also forces acting 
on line, so called line forces, can exist in nature and are typically surface tension but is 
not considered here. The inertia forces come from the acceleration of a fluid element, or 
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APPENDIX A. FLUID FLOW EQUATIONS 



advection, and are expressed in terms of the material derivative 



Dili 

~Dt 



dui duj 
dt 1 dx j 



The derivative describes the change in velocity as the element moves in space. 



(A. 3) 



A.2 Conservation of Momentum 



By applying Newton’s law of motion to an infinitesimal fluid element the equation of 
motion relating the inertia to the net force at a point for any continuum becomes 



Duj dTjj 
P Dt dxj 



(A.4) 



which is the Cauchy’s equation of motion. The constitutive equation for a particular fluid 
relates the shear stress to the strain rate. This relation is embedded in the stress tensor and 
could take a complex form. Many fluids may be approximated by a Newtonian model. 
The Newtonian model uses the assumption that the shear stress is linearly proportional to 
the strain rate. When the rate of shear strain is zero, the shear stress is zero. 



A.3 Navier-Stokes Equations 



The Navier-Stokes equations are obtained by substituting a constitutive relation into Eq. (A.4) 
and by assuming a Newtonian fluid becomes 



P 



Du; 

~Dt 



dp_ + d_ 
dx/ dxj 




|l(V.u)5y 



(A. 5) 



which is a general form of the Navier-Stokes equations where T| and p could be functions 
of for example pressure and temperature. 8 is the Kronecker delta and the strain rate tensor 
e is defined as 



1 / diii du i 

2 V 3* i dxj 



(A. 6) 



If the Navier-Stokes equations are written in a three dimensional coordinate system they 
are three in numbers but include four unknowns, u and p. To be able to solve the problem, 
a complement to the conservation of momentum the principle of conservation of mass can 
be employed, called the continuity equation: 



dp 

dt 






(A.7) 



If the variation in viscosity is small the Navier-Stokes equations can be reduced to 



A .4. STOKES EQUATIONS 
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If the variation in density is small, Vu = 0 and the Navier-Stokes and the continuity equa- 
tions reduce to 

Dm _ 

p— = -Vp + TlV-U 

V u = 0. 

If the no variation in time occurs the Navier-Stokes equations can be further reduced: 

p(u- V)u = -V/7 + r|V 2 u. (A.9) 

If the viscous effects (the rightmost term) are negligible the equations can be reduced to 
the Euler equations: 

p(u • V)u = — Vp. (A.10) 



A.4 Stokes Equations 

In lubrication problems however the fluid shear stress is normally high due to the small 
film thickness and cannot be neglected. On the other hand, since the shear stress often 
is high in hydrodynamic lubrication the fluid inertia or advection often becomes small in 
comparison. Hence, a more realistic simplification would be to neglect the advection: 

— V/7 + r |V 2 u = 0 (A. 11) 

which are known as the Stokes equations. These are more convenient equations to handle 
since the non-linear advective term is dropped. These equations are well suited for creep- 
ing flows with large viscosity and small velocities. In Paper A, Navier-Stokes Eq. (A.9) 
and Stokes Eq. (A. 11) are compared in a hydrodynamic lubrication application. 



A.5 Reynolds Equation 



The Navier-Stokes equations can be reduced further by assuming that the length scales 
L y /L x <^i 1 see Fig. 2.1. Consider the Stokes equations (A. 1 1) written in two dimensions: 



dp 

-fc + " 

dp 



iijj, fy\ 
V dx 2 dx~ ) 

A+A 

W dy-J 



0 . 



We will consider a scale analysis which could equally well be applied in three dimen- 
sions but for simplicity reason we consider a two-dimensional case here. Introduce the 
dimensionless parameters: 



t 2„ 

x y u v L y p 

x = — ; y = — ; u = — ; v = — ; p = — : 

L x Ly li {] Vq L x T\qUq 



(A. 12) 
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Inserting these dimensionless parameters into Eq. (A. 12), dividing by (r|o«o) and multi- 
plying by L,] gives: 



dx \L X J \dx 2 u 0 dx" J 

L x dp _ (d 2 u v 0 3 2 v\ 

d + B —v H t = 0. 

Ay dy \ dy «0 dy~ J 



The ratio vo /mo has the same order of magnitude as L y j L x . Neglecting these terms and 
assuming that that pressure variation in the y-direction can be neglected the equation be- 
comes: 



dp 

aT 1 ^ 



. d 2 u 

% 2 



(A. 13) 



and the Navier-Stokes equations have been reduced by applying the thin film approxima- 
tion. This equation is in dimensionless form. However, let us express the equation with 
dimensional quantities and integration twice with respect to y yields 



1 dp o 

U= 2^Tx r+Ay + B - 



(A. 14) 



A and B in this case are either constants or functions of x and can be solved by applying a 
boundary condition for u: 



u = Uo at y = 0 and 

u = Uh at y = h. (A. 15) 



where h = h(x). Substituting (A. 14) into (A. 15) gives 

1 . , h dp 

B = U 0 and A=-(U h -Uo)-—rf. 

h 2r| dx 

The Velocity becomes 

« = ~y h ) +l(Uh - u 0 )+u 0 . 

Revisiting the continuity equation for incompressible conditions: 

du dv 

a^ + a^ = 

and integration with respect to y from 0 to h(x) yields: 



(A. 16) 



(A. 17) 



(A. 18) 



rK x ) du 

-(Vo -V h )=x-j ^ dy (A. 19) 

By substituting (A. 17) into (A. 19) interchanging the differentiation and integration through 
the use of applying Leibnitz’s rule we get 



-(Vo -VO = - 



-(Vo -VO = 



1 dp /■*(*) 
2r| dx 
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dx 
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The left-hand term —(Vo — V/,) = dh/dt. In this case dh/dt = 0 because of steady state 
conditions and the equation can be written: 



This is the incompressible Reynolds equation describing fluid film flows where the scales 
Lyl L x 1 and pressure variations across the film is neglected. The equation is here one 

dimensional with the y-componcnt transformed to the film thickness equation h depending 
only on x. If the coordinate z was not dropped the film thickness would be a function of 
both x and z. In either case, all the Navier-stokes equations and the continuity equation are 
reduced to one equation with pressure as the dependent variable. If U], zero and u = Uq 
the Reynolds equation is reduced to 



This equation is used as a basis and modified in different ways in Paper B and C. 




(A. 20) 




(A. 21) 



In a similar form the compressible Reynolds equation can be written 




(A. 22) 
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Abstract 

Results of a numerical study of the influence of micro-patterned surfaces in hydrodynamic 
lubrication of two parallel walls are reported. Two types of parametrized grooves with the 
same order of depth as the film thickness are used on one stationary wall. The other wall 
is smooth and is sliding with a specified tangential velocity. Isothermal incompressible 
two dimensional full film fluid flow mechanics is solved using a Computational Fluid Dy- 
namics method. It is shown that, by introducing a micro-pattern on one of two parallel 
walls, a net pressure rise in the fluid domain is achieved. This produces a load carrying 
capacity on the walls which is mainly contributed by fluid inertia. The load carrying ca- 
pacity increases with Reynolds number. The load carrying capacity is reported to increase 
with groove width and depth. However, at a certain depth a vortex appears in the groove 
and near this value the maximum load carrying capacity is achieved. It is shown that the 
friction force decreases with deeper and wider grooves. Among all geometries studied, 
optimum geometry shapes in terms of hydrodynamic performance are reported. 
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A.l Introduction 

The purpose of lubrication is to separate surfaces in relative motion in order to reduce 
friction and wear. The separation and load carrying capacity are achieved by generating 
a pressure in the fluid film between the surfaces. The most significant pressure build- 
up in hydrodynamic lubrication is achieved when a converging gap is allowed to form 
the surfaces. A relative velocity between the surfaces will start driving fluid into the gap 
leading to a pressure build-up. 

In order to improve performance and efficiency and to reduce wear and the risk of fail- 
ure in hydrodynamic lubrication, it is important to study various effects that can contribute 
to the pressure build-up. 

The subject for this work is to investigate the influence of surface gesturing in hy- 
drodynamic lubrication of two parallel surfaces. Experimental research in this area, such 
as reported in [1, 2] shows that significant improvements in hydrodynamic journal bear- 
ing performance can be achieved by surface gesturing. Textured reciprocating surfaces are 
studied experimentally in [3] where a clear reduction in friction force is found. Glavatskih, 
et al. [4] investigated hydrodynamic performance of thrust bearings and conclude that 
higher film thickness and lower power loss are obtained through the use of micro-patterned 
pads. 

All mechanisms contributing to bearing efficiency using surface gesturing are far from 
known. Therefore theoretical analysis is important to be able to study individual effects. 
Numerical studies are presented in [5, 6] where Reynolds equation is used to simulate the 
effects of spherical grooves applied on a stationary surface in hydrodynamic lubrication. 

One effect contributing to the load carrying capacity could be local cavitation at the 
diverging part of surface grooves. Local cavitation could contribute to the pressure build- 
up in a bearing. Such analysis is done by Brizmer et al. [7] who study parallel thrust 
bearings with spherical grooves using Reynolds equation with a cavitation model. 

In the diverging part of the groove cavitation may occur leading to a two-phase flow. 
This two-phase flow will then enter the convergent part of the groove. To achieve a net 
pressure build-up this flow must be transformed to a single -phase liquid flow. Because of 
the groove and film geometry, this can only be achieved if an additional amount of liquid 
is supplied from outside the groove. In practical cases, it is unclear how this liquid supply 
to the dimples is achieved, especially when side leakage is considered. That is why there 
is a possibility that cavitation has a limited effect. Thus, other effects that may lead to a 
pressure build-up must be studied. 

Arghir, et al. [8] report a Computational Fluid Dynamics (CFD) analysis of fluid flow 
between parallel surfaces with relative speeds. One of the surfaces contains a groove which 
is represented by three different geometries. The main conclusion of this study is that there 
is a net pressure build-up caused by the combined effect of the presence of a groove and 
of the increasing Reynolds number. 

Although there are many physical effects contributing to hydrodynamic efficiency, it is 
important to study them individually to be able to distinguish between them. For different 
applications and running conditions, different physical effects will be present, contributing 
to the hydrodynamic performance in various ways. In this paper the fluid mechanical 
effects will be isolated and analyzed. The importance of inertia for certain types of flow 
conditions and micro-patterns will be discussed further in this work. The hydrodynamic 
performance in terms of friction force and load carrying capacity and how it depends on 
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groove geometry and flow conditions will be studied by using CFD. 



A.2 CFD-Model 



In this work the Navier-Stokes equations are solved using a commercial CFD code. The 
equations are applied without body force, with constant viscosity and density. The equa- 
tions are steady and solved in the x- and y-direction only. With these properties the Navier- 
Stokes and the continuity equations can be written, respectively: 

p(u-V)u = -Vp + r|V 2 u 
V u = 0. 



The equations are solved in dimensional form, however, all results are represented in 
non-dimensional form in such way identical result are achieved for the same Reynolds 
number, Re. Hence, the dimensional values do not have any meaning, only the value of 
Re. Information about the non-dimensionalization of the quantities can be found in [9]. 

In this work, two types of geometries for the fluid domain have been studied. The first 
geometry is defined as the cylindrical geometry and is seen in Figure A. 1 . The groove is 



y 



0 - 




Figure A.1: The fluid domain of the cylindrical geometry. Included are the geometrical 
parameters and boundary conditions. 

described by a circular arc of depth d and width vv. The constant quantity L x = \e — 3 m is 
the length in the x-direction of the domain and the constant quantity L y = 3e — 5 m is the 
film thickness of the fluid domain. In the same figure, the boundary conditions for the fluid 
domain are shown. The fluid is experiencing no-slip conditions at the walls. The upper 
wall is moving with a velocity u whilst the lower wall is stationary. The edges at high and 
low x-coordinates are bound by a periodic boundary. A symmetry condition is applied in 
the z-plane. 

The second geometry is defined as the splined geometry and is seen in Figure A.2. The 
geometrical parameters for this case are defined as for the cylindrical geometry except that 
the groove is described as a splined curve. Another parameter, xd , is introduced to describe 
the x-displacement of the low y point of the groove. The same boundary conditions are 
applied for the splined as for the cylindrical geometry. 

In all the equations, terms are discretized in space using second-order centered differ- 
encing apart from the advection. For these terms, a hybrid differencing scheme is used 
which is close to first-order accuracy [10]. This scheme is slightly better than upwind 
differencing because second-order central differencing will be used across streams and in 
regions of low flow where Pe < 2. If Pe > 2, upwind discretization will be applied. These 
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Figure A. 2: The fluid domain of the the splined geometry. Included are the geometrical 
parameters and boundary conditions. 



discretizations are based on a conservative finite-volume method [11]. For information of 
the computational grids and error analysis, see Sahlin [9]. 

A parametrized fluid model is used where certain parameters are being changed in a 
systematic order. All the parameters that are being varied, entirely independently from 
one another, are listed in Table A.l. Results from some complementary simulations with 



Table A.l : Individually varied parameters. 



'0.15 
0.20 
0.25 
0.30 
* 0.35 
0.40 
0.45 
0.50 



'0.25 




0.50 


f-0.3 


0.75 xd+ < 


0 Re < 


1.00 


+0.3 


1.25 





40 

80 

120 

160 



values different than those in the table will occasionally be shown. 

A.3 Results and Discussion 

Results are presented and discussed in terms of advection, fluid film pressure, streamlines 
and forces acting on the walls. It should be emphasized that the fluid flow effects are com- 
pletely isolated from any other effects that might occur in a real thin film fluid flow. The 
interest in this paper is thus to analyze the fluid mechanical effects in hydrodynamic per- 
formance. However, other effects such as cavitation and thermodynamics would certainly 
influence the hydrodynamics. 

A.3.1 Advective Influence 

The advective (convective) contribution in the fluid flow can be studied by comparing the 
Navier-Stokes solution with the Stokes where the advective terms (inertia) are truncated 
from the Navier-Stokes equations. In Fig. A.3 the pressure distribution on the upper wall 
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(a) Re = 10 




(b) Re = 160 



Figure A.3: Comparison between Navier-Stokes and Stokes solutions for pressure distri- 
bution on the upper smooth wall for the cylindrical geometry. w + = 0.2 and cl + = 0.25 for 
both plots. The vertical dotted lines represent the grove edges. 



has been plotted as a function of the spatial dimension x + for Navier-Stokes and Stokes 
solutions respectively. In this figure the cylindrical geometry is used where Re is set to 
10 in A. 3(a) and 160 in A. 3(b). With the value of Re in Fig. A. 3(a), the Navier-Stokes 
and Stokes solutions are much alike. A small displacement to the left of the Navier-Stokes 
curve is seen next to the portion of the groove. When Re is increased in Fig. A. 3(b), great 
differences between the two solutions arise. From Fig. A. 3(b), the Stokes curve seems to 
be entirely antisymmetric whereas the Navier-Stokes curve shows a net pressure rise in the 
domain from the zero boundary. 

Following the pressure curves in Fig. A. 3(b) starting from low x + , an almost linear 
pressure drop is encountered for both curves. Slightly before arriving at the groove edge, 
located at x + = —0.2, pressure increases rapidly for the Navier-Stokes solution producing 
a net pressure build-up in the fluid film. For the Stokes solution, the pressure continues 
its descent until it passes the groove edge where the curve changes to an almost linear 
pressure rise. The solutions from the Stokes case produces close to an antisymmetry-axis 
at the line x + — 0. 

The area between the two curves in Fig. A. 3(a) as well as in A. 3(b) is equal to the 
contribution of inertia in the flow only, because the advective term is completely dropped 
in the Stokes equations. This confirms findings of a recent study made by Arghir et al. [8]. 

A quantitative measure of the difference between the Navier-Stokes and Stokes So- 
lutions can be seen in Fig. A. 4, which shows the dimensionless normal force F+ as a 
function of w + for Navier-Stokes and Stokes solutions. As an example, the cylindrical 
geometry with the lowest simulated groove depth, d + = 0.25, is chosen. This is in order 
to provide as low an influence of inertia as possible. In Fig. A. 4(a) Re = 10. 
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(a) Re = 10 



(b) Re = 160 



Figure A.4: F + on a 10-base logarithmic scale as a function of w + for Navier-Stokes and 
Stokes solution respectively, d + = 0.25. 



It is seen that there is indeed a positive net force in the direction normal to the wall for 
the Stokes solution, evolving from diffusion. 

However, the force for the Navier-Stokes solution is one order of magnitude greater. In 
Fig. A. 4(b) the value of Re is 16 times greater and the force for the Navier-Stokes solution 
is two orders of magnitude greater than the force from the Stokes solution. 

This shows the importance of the advective terms in the Navier-Stokes equations for 
this type of geometries and flow conditions. It implies that further simplifications to the 
Navier-Stokes equations would not be appropriate here. 

A.3.2 Pressure Distribution 

Figures of pressure distributions are shown with values at the periodic boundary set to 
zero only because it gives a convenient representation. The incompressible Navier-Stokes 
equations are solved, hence, the absolute pressure has no meaning in this study and the 
pressure distribution across the domain could be shifted any amount. Only the gradient of 
the pressure is affecting the flow and has meaning in this case. 

Comparisons of pressure distribution for the cylindrical geometry are made in Fig. A. 5 
where p + is plotted as a function of the non-dimensional spatial dimension x 1 for different 
values of Re. w + = 0.3 is constant for all sub-figures but d ' is varied between the sub- 
figures. All pressure distributions show a linear pressure drop at the beginning and at the 
end of the domain. The length in the x-direction where the pressure is deviating from this 
linearity at the mid part of the domain will be called the affected elongation. If Fig. A. 5(a) 
where d + = 0.25 is considered, it is clear that the pressure build-up in the groove increases 
with the value of Re, and the differences are substantial. The pressure curves are con- 
vexly shaped across the entire affected elongation, i.e., the curvature d 2 p/dx 2 is negative 
throughout the pressure rise for all values of Re. With Re = 160, the maximum dimension- 
less pressure is about 9. In Fig. A. 5(b) the groove depth is increased to d + = 0.75 and it 
is seen that the maximum non-dimensional pressure has increased to about 12. However, 
the curves do not show as extreme convexity as in the previous figure. In other words. 
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(a) d + = 0.25 




(b) d+ = 0.75 




(c) d + = 1 .25 



Figure A.5: Upper wall pressure distribution for the cylindrical geometry with w + = 0.30. 
The vertical dotted lines represent the grove edges. 



the curves are flattened out in the affected elongation. In Fig. A. 5(c) the greatest groove 
depth d + = 1.25 is used. Here the pressure maximum is decreased compared to the case 
of d + = 0.75 and the difference in solutions between the Reynolds numbers are decreased. 
At the right half of the affected elongation the pressure distributions have concave shapes. 

Some preliminary conclusions can now be made regarding how the pressure distribu- 
tions are affected, which will later be proved with further results. By studying Fig. A.5 it 
seems that the affected elongation depends on the groove width w + . On the other hand, 
the value of d + seems to affect the elevation and form of the pressure build-up curve. In 
terms of load carrying capacity, a convexity in the affected elongation which gives the 
maximum possible area beneath the curve is preferred. A flat pressure build-up with a 
higher maximum pressure will not lead to a higher load carrying capacity. 

In Fig. A. 6 the pressure distributions are plotted for the splined geometry with w + = 
0.5 and d + = 0.25. Fig. A. 6(a) is for a negative displacement of the low y-point with xd + = 
—0.3, Fig. A. 6(b) for xd + = 0 and Fig. A. 6(c) for a positive displacement with xd + = 
+0.3. As the affected elongation is nearly the same for all sub-figures, it can be confirmed 
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(a) xd + = -0.3 




(b) xd + = 0 




(c) xd + = +0.3 



Figure A.6: Upper wall pressure distribution for the splined geometry used with d + = 0.25 
and w + = 0.50. The vertical dotted lines represent the grove edges. 



that it strongly depends on the value of w + . It seems that the pressure distributions in 
Fig. A. 6(a) provide highest load carrying capacity. When the groove is symmetric, in 
Fig. A. 6(b), the rising-pressure zone flattens out slightly which may lead to lower load 
carrying capacity. In Fig. A. 6(c) where xd + = 0.3, the rising-pressure zone flattens out 
further but a slightly increased pressure maximum is noted. If the assumptions about the 
load carrying capacities are true, the magnitude of the pressure maximum in the domain 
is not of primary importance in the sense of hydrodynamic performance. More important 
would be to achieve a shape of the pressure distribution which is as much convex at the 
affected elongation as possible. 



A.3.3 Streamlines 

Streamlines can give a qualitative idea about how the flow is ordered, if any vortexes or 
irregularities of the flow exist and, in that case, where they develop and what influences 
they have on the hydrodynamic performance. Streamlines are expressed in terms of the 
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scalar stream function i|/ and along a streamline \|/ = constant. The velocities can be found 
in terms of the derivatives of \|/ as u = ()\\i /dy and v = — 3\| f/dx. 

In Fig. A. 7, streamlines for the cylindrical geometry are plotted, all with Re = 40 




(d) d+ = 1.25. 



Figure A.7: Streamlines in the cylindrical geometry for different values of d + are plotted 
where Re = 40 and w + = 0.20. 

and groove width w + = 0.20 but different groove depth d ' . In A. 7(a) the streamlines 
are completely ordered and smooth, but as d + increases to a value between 0.5 and 0.75, 
a vortex is developed in the groove (see Sahlin [9] for a comparison with isobar plots 
for the same conditions). The vortex development can be related to a driven cavity flow 
phenomenon extensively studied by Shankar and Deshpande [12]. The size of the vortex 
is increased with increased d + and at d + = 1 .25 the vortex is about the same size in the 
y-direction as the groove depth value. This shows that the vorticity is clearly dependent on 
groove depth. 

Figure A. 8 shows the effects of different values of Re. The cylindrical geometry is 
used with w + = 0.15 and d + = 0.25. As the value of Re increases, the size of the vortex 
increases. There is no vortex at a sufficiently low value of Re. 

In Fig. A. 9, streamlines are plotted for the cylindrical geometry with constant groove 
depth d + but different w + . It can be seen that the vortex decreases in size with increasing 
groove width w + and constant Re. At a sufficiently large value of w + there is no vortex 
and no back-flow exists in the groove. 

The vortex development obviously depends on the value of Re and how fast the ge- 
ometry is changing, i.e. the value of dy/dx. These are properties which highly affect 
the advective influence whereby a connection is likely to exist between advection and the 
amount of vorticity. Let us therefore compare the full Navier-Stokes equations with the 
Stokes in terms of streamlines. Such a comparison is made in Fig. A. 10. By comparing 
A. 10(a) and A. 10(b) it is obvious that the inclusion of advection leads to an earlier devel- 
opment of the vortex. However, when d + is doubled in A. 10(c) and A. 10(d), it can be seen 
that a vortex develops even without the advective influence and is entirely centered in the 
x-direction. With the advective influence the vortex center is moved to the right. Hence, 
the streamlines for the Stokes solution are completely symmetric in the line x + = 0. As 
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(d) Re = 160 



Figure A.8: Streamlines in the cylindrical geometry for different values of Re are plotted 
where d + = 0.25 and w + = 0.15. 




(a) w+ =0.15. 




(b) w+ = 0.25. 



(c) w+ = 0.35. 



(d) w + = 0.45. 



Figure A. 9: Streamlines in cylindrical geometry for different values of w + . Re = 40 and 
d + = 0.75 in all figures. 



shown in the pressure plots, Fig. A. 3, the pressure distributions for the Stokes solutions 
are close to antisymmetric at the same position. 

A.3.4 Forces 

Wall forces consist of both inertial and diffusive forces. These are attained directly from 
the CFD-solver. 

It was shown in Fig. A. 5 and A. 6 that the pressure build-up was monotonically in- 
creased by increasing Reynolds number. This increase is further substantiated by Fig. A.l 1 
where F. + is plotted as a function of w 1 . It is shown that the normal force also monoton- 
ically increases by increasing groove width w + . A quantitative comparison between the 
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(a) w f = 0.15, d + = 0.5. Advection. (b) w f — 0.15, d + = 0.5. No advection. 




(c) w + = 0.15, d + = 1. Advection. 




(d) w + = 0.15, d + = 1. No advection. 



Figure A. 10: Streamlines for Navier-Stokes solution left and Stokes solution right. The 
cylindrical geometry is used and Re = 160 for all. 




(a) Cylindrical, d + = 0.25 




(b) Splined, d + = 0.25 



Figure A.1 1 : F+ as a function of w + for various values of Re for the cylindrical geometry 
(top) and the splined with xd + = 0 (bottom). 



cylindrical geometry and the splined geometry with a symmetrical groove displacement 
xd = 0 shows that higher normal forces are achieved with the cylindrical geometry for 
d+ =0.25. 

Figure A. 12(a) shows how the normal force is affected by different values of xd ' for 
the splined geometry. As w + increases, the differences of /y between the three values 
of xd + also increase. It is clear that a negative value of xd 1 contributes positively to F v + , 
especially when w + is large. By comparing this curve with Fig. A. 11(a) for Re = 40 it 
is noted that the splined geometry with xd + = —0.3 produces higher values of F v + that 
the cylindrical. It should also be noted that F v + , for the special case of xd + = 0.3 in 
Fig. A. 12(a), does not monotonically increase as w + increases. This is an exceptional 
behavior compared to all other geometries. If F + is instead plotted as a function of d + , 
the differences are larger when d + is small and decrease when d + increases. This can be 
seen in Fig. A. 12(b). It can also be seen that an optimum value of d + exists on each curve. 
There is a tendency that the optimum is displaced towards greater values of d + as xd + 
increases. 
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w + d + 

(a) Re = 40, d + = 0.25 (b) Re = 40, w + = 0.15 



Figure A. 12: F + for various values of xd + for the splined geometry.. 

In Fig. A. 13 the non-dimensional friction force is plotted against d + for the cylindri- 
cal and splined geometry respectively. The cylindrical geometry produces lower friction 





Figure A. 13: F+ as a function of d + for the values of w + = {*0.15, >0.2, +0.25, DO. 3, 
x0.35, 00.4, V0.45, A 0.5}. Solid line represents Re = 40 and dashed line Re = 160. 

force than the splined geometry for all cases. This is partly explained by the fact that the 
cylindrical groove inherits a greater area than the splined, using the same width and depth, 
which in turn gives rise to a lower friction force. It is noticed that F x differs between Re 
and the difference is increasing as d + increases. The difference in Fj between Re does 
not, however, change much by different values of w + . F x decreases with increasing val- 
ues of w + and d + . However, Fj shows a different behavior for the cylindrical geometry 
with high Re and low w + . The smooth surface case, i.e. d + = 0, corresponds to a value of 
F x = 3 for a value of Re = 40. This shows that the textured surface will give rise to lower 
friction. 



A.3. RESULTS AND DISCUSSION 



37 



It was shown in Fig. A. 1 1 that the cylindrical geometry produces greater normal force 
which is true for some cases. It was concluded in Fig. A. 12 that the splined geometry 
is most effective, in terms of normal force, with a negative value of xd + . Figure A. 13 
showed that the cylindrical geometry produces lower friction force compared to the splined 
geometry. The differences in friction force for different values of xd + is negligible. This 
concludes that, in most cases, the most beneficial geometries in terms of hydrodynamic 
lubrication are the cylindrical geometry and the splined geometry with xd + = —0.3. These 
two geometries are compared in Fig. A. 14 and A. 15, where the normal force is plotted as 
a function of d + for different values of w 1 . At the lowest Reynolds number in Fig. A. 14, 





d + d + 

(a) Cylindrical. Re = 40 (b) Splined, xd + = —0.3, Re = 40 

Figure A. 14: F + as a function of d + for Re = 40 and the values of w + = {*0.15, >0.2, 
+0.25, DO. 3, x0.35, 00.4, V0.45, A 0.5}. The cylindrical geometry (top) and the splined 
(bottom). 




Figure A.15: F + as a function of d + for Re = 160 and the values of w + = {*0.15, >0.2, 
+0.25, DO. 3, x0.35, 00.4, V0.45, A 0.5}. The cylindrical geometry (top) and the splined 
(bottom). 

the splined geometry shows a higher load carrying capacity. At Re = 160 in Fig A.15 
however, the properties are reversed and the cylindrical geometry shows higher normal 
forces. As stated earlier, there exists an optimal value of d + in terms of normal force. 
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which is also the case in this figure. At the lowest value of w + the maximum is located 
at d + 0.5 in all cases represented. As w + increases, the maximum is displaced towards 
greater groove depths. This is also the case when Re is decreased. However, the optimum 
could in general be fit in the range 0.5 < d + < 0.75 for most cases. No direct relation 
between the ratio d/w and forces can be found because of the strong individual effects of 
d + and w + respectively. 

The conclusion that a certain maximum groove depth exists where an increase does 
not lead to further increase in load carrying capacity agrees with experiments made by 
Snegovskii and Bulyuk [1]. 

A comparison could be made between the second curve from the bottom representing 
w + = 0.20 in Fig. A. 14(a) and the streamline plots in Fig. A. 7(a) and Fig. A. 7(b) which 
represent the same geometry and Re. It is seen that at d + « 0.5 where the maximum 
load carrying capacity is located, a vortex is beginning to form. When Re is increased, 
the vortex and the Pj -maximum appears at lower values of d ' which is evident from 
Fig. A. 8 and A. 14 respectively. Making a comparison between Fig. A. 9 and Fig. A. 14 and 
A. 15 reveals the relation between vortex formation and fy -maximum in terms of different 
values of w + . 

This discussion shows that the normal force maximum appears close to the beginning 
of a vortex formation in the groove. This is true for most cases studied. Hence, the vortex 
affects the normal force in a destructive manner. 



A.4 Conclusions 

The fluid mechanics of hydrodynamic lubrication between parallel surfaces is studied for 
isothermal, incompressible and steady 2D conditions. The conclusions are summarized as 
follows: 

1 . An introduction of a micro groove on one of the surfaces affects the flow and pres- 
sure pattern. This gives a net pressure build-up and a load carrying capacity of the 
film. Negligible hydrodynamic effects are achieved when the advective terms are 
excluded from the Navier-Stokes equations. 

2. Load carrying capacity increases with Reynolds number and groove width w + . A 
general optimum is achieved when 0.5 <d + < 0.75 for all geometries and groove 
widths. The load carrying capacity can, to some extent, be related to the amount of 
circulation in the groove. A vortex appears at a certain value of groove depth d + . 
Near this value the maximum load carrying capacity is achieved. 

3. The friction force decreases with increasing values of groove depth d + and width 



4. The best hydrodynamic performance is achieved with the cylindrical geometry and 
the splined with a negative groove displacement xd + . 
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Abstract 

A general cavitation algorithm is presented that accommodates for an arbitrary density- 
pressure relation. It is now possible to model the compressibility of the lubricant in such 
a way that the density-pressure relation is realistic through out the contact. The algorithm 
preserves mass continuity for cavitation caused by bearing geometry as well as surface 
topography. It is a commonly accepted physical assumption that the contribution of the 
pressure driven flow can be considered negligibly small in the cavitated region. This phe- 
nomenon is adopted in the present algorithm, which is similar to that of Elrod, and is 
modeled by a switch function that terminates the pressure gradient at the regions of cavita- 
tion. Results with this algorithm for different density-pressure relations are presented and 
discussed. Effects of inlet conditions, such as surface roughness and starvation, on load 
carrying capacity of the contact are analyzed. 
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B.l Introduction 



Cavitation usually occurs in the regions of diverging contact gap and implies sub-ambient 
pressures. These low pressures lead to a transformation of the liquid into a gas-liquid mix- 
ture. Different types of cavitation models have been proposed over the years to replicate 
this behavior in theoretical simulations. Some of the algorithms preserve lubricant mass 
continuity whilst others not. By simulating the behavior of rough surfaces, several cav- 
itation regions might appear inside the contact and preserving mass continuity becomes 
crucial. 

Jakobsson and Floberg [13] developed a mass preserving cavitation theory. They as- 
sumed that the pressure is constant in the cavitation region which means that the pressure 
gradient is zero. They derived a set of conditions to locate the cavitation boundaries. 
Later, Floberg extended their theory of cavitation and implemented it for a number of 
bearing types [14, 15, 16, 17]. Elrod [18] and Elrod and Adams [19] developed a cavita- 
tion algorithm where a single equation is used throughout the lubrication region without 
the need for explicit equations to locate the cavitation boundaries. They used a switch 
function to terminate the pressure gradient in the region of cavitation. Vijayaraghavan 
and Keith Jr [20] developed an algorithm in a similar way as Elrod and introduced more 
rigorous derivation which was further developed in [21]. This algorithm is derived under 
the assumption that the bulk modulus is constant. For real lubricants the bulk modulus 
varies with pressure. For a wide pressure range and under certain conditions it becomes 
important to use a realistic compressibility model for the lubricant. However, treating the 
bulk modulus as a constant could produce good results for a narrow pressure range. 

The cavitation algorithm presented here uses Reynolds equation and is based on the 
assumption that the pressure gradient is zero in the cavitation region and that only Couette 
flow is present there. It inherits a similar switch function as the Elrod algorithm to modify 
Reynolds equation at the locations of cavitation. With this algorithm, an arbitrary density- 
pressure relation of the lubricant can be applied. Thus, measured compressibility factors 
from real lubricants can be used to achieve more realistic results. 



The governing equation that constitutes the basis for the current algorithm is the Reynolds 
equation. A simple steady state form of this equation can be written as 



This is a mass continuity equation with the pressure p as the solution variable, describing 
Couette and Poiseuille flows, the first and second term respectively. It is non-linear but 
becomes linear if h, p and q are independent of p. The Reynolds equation is well suited 
for thin film flow, where h, and rate of change in h, is small with respect to the other spatial 
coordinates. 

The density of all liquid lubricants depends on pressure among other variables. The 
density-pressure relation changes abruptly when the lubricant undergoes cavitation, re- 
sulting in a mixture of liquid and gas. At this condition the pressure remains close to 
constant whilst the density decreases rapidly. In order to mimic the behavior of cavitation, 
the Reynolds equation needs to be modified. It is convenient to use the density as the 
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(B.l) 
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dependent variable instead of pressure. We start by defining the non-dimensional density 
as 

e = — (B.2) 

pc 

where p c is the density at ambient pressure p c . Since p = /? ( 0 ) it is possible to apply the 
chain rule on the pressure gradient, i.e., 

V/7= ^ V9 ' (B3) 



Combining Eqs. (B.l), (B.2) and (B.3) we get: 



m 3(0/0 _ v / 9 /0 dp \ 

2 dx \ 12r| d0 J 



= 0 



(B.4) 



where the non-dimensional density 8 is now the dependent variable and the only pressure 
influence in the equation is through the pressure derivative d/; /d0. In the cavitation region 
the pressure is assumed to be constant, p = p c , [13]. Hence, the pressure gradient is zero: 



Vp = 0, 



(B.5) 



and no Poiseuille flow contribution exists in this region. However, the mass will be pre- 
served and transported only by means of Couette flow. By definition, 9 = 1 at the cavitation 
boundary, thus to accomplish the termination of the pressure gradient we utilize the unit 
step function 



j 1 , 9 > 1 ; Full film region 

\ 0, 9 < 1; Cavitation region. 



Applying the above step function to the expression for the pressure gradient in Eq. (B.3) 
yields the Reynolds equation in the form 



u 3(9 h) 
2 dx 



( — s^ve 

\12r| 00 



= 0. 



(B.7) 



Here, a single equation which preserves lubricant mass describes the flow both in the 
full film and in the cavitation region respectively. This is a more general form of al- 
ready existing cavitation algorithms based on the Reynolds equation in that an arbitrary 
density-pressure relation could be used. Thus, more realistic density-pressure relations 
that accurately reflect the behavior of real lubricants can easily be applied. 



B.3 Compressibility 

The compressibility of the lubricant can be treated in several ways. It can be expressed 
through the bulk modulus defined as 

e= e K (B - 8) 

If [3 is considered constant, integration ones gives the the expression for pressure and 
density: 
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This expression for constant bulk modulus is sometimes used to model the compressibility 
of lubricants. However, in real lubricants (3 varies with density and thus the expression is 
only valid in a limited pressure range. Another is the Dowson & Higginson expression 
[ 22 ]: 



Q +C 2 (p — p c ) 
Cl + (p- Pc) 



(B.10) 



where the coefficients C\ = 0.59* 10 9 and C 2 = 1.34 was fitted to measured compress- 
ibility data of a mineral oil. Since the Reynolds equation on the form presented here can 
obey an arbitrary compressibility relation for the lubricant, it would be reasonable to use 
some other measured data as well as an input. Hence, as a test case for the current cavita- 
tion algorithm, the compressibility for a mineral oil will be used where data is taken from 
Tuomas and Isaksson [23]. They carried out measurements in a high pressure chamber 
with pressures ranging from ambient up to about 3 GPa. The coefficients in the Dowson 
& Higginson expression where fitted to the measured data. This fit can be seen in Fig. B. 1 
together with the measured data and the original Dowson & Higginson expression. The 




Figure B.1 : Compressibility models for the measured mineral oil data. The original Dowson 
& Higginson expression with C\ = 0.59e9 and C 2 = 1 .34, a Dowson & Higginson expression 
fitted to the measured data with C\ = 2.22e9 and C 2 = 1 .66, an exponential expression with 
(3 = 3.3 GPa fitted at 9 = 1 for the fitted D & H expression. 

relation in Eq. (B.9) is also plotted, where (3 = 3.3 GPa which is the bulk modulus taken 
for the curve fit at 0 = 1. The best correlation between the expressions can be observed in 
the proximity of ambient pressure but deviates exponentially with pressure rise. The same 
notation for the compressibility relations in Fig. B.l will be used in the figures ahead. The 
acronym MNL is used as a label to show that the compressibility relation origins from 
the measured mineral oil data. The cavitation pressure p c = 10 5 Pa will be used in the 
numerical examples. Cavitation and ambient pressure will be treated synonymously. 



B.4 Numerical Examples 

The Reynolds equation is discretized by using a central differencing on the Poiseuille term 
and a second order backward differencing on the Couette term. Because of the non-linear 
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nature of the problem, a fixed point iteration process on the form 0„ = /(0„_ 1 ,...) for 
n = 0,1,... is used to solve the discrete system. The switch function g is updated each 
iteration. The iterative convergence of the solution is always ensured to be a value below 
1 0 s and the grid is ensured to be fine enough to produce a grid independent solution. 

In order to compare solutions for different compressibility models etc., the following 
integral residual of the current solution \|/ to a solution \|/* is defined as: 





rL x rL x 




I Lx 


c* = 

‘ J, \j f 


/ V(/*dA - / \|/d.v 

JO JO 


/ 


/ \|/ dx 

Jo 



which is the residual in load carrying capacity if the solution is pressure. In the numeri- 
cal examples, only one dimensional problems will be considered where L x is the domain 
length. This means that the geometries are treated as infinitely wide with no side leakage. 

Dirichlet boundary conditions are used at the inlet (w) and the outlet (e) of the domain 
by specifying a constant pressure p or non-dimensional density 0. If 0 < 1 at the boundary, 
the domain will be termed starved, likewise if 0 = 1 the domain is flooded and if 0 > 1 the 
domain is pressurized. As a first comparative test case, a parabolic geometry modeled by 

H = H 0 +A ^~ Lx ^ (B.12) 

(Lx/2) 2 

with the same physical dimensions and operating conditions as in Vijayaraghavan and 
Keith Jr [20], the geometry is seen in Fig. B.2. In Fig. B.3, a pressurized inlet is used 




°0 “ *■ 0.0762 



Figure B.2: A single parabolic slider geometry with Hq = 25 pm and A = Hq. 



and the pressure solutions are shown from the Dowson & Higginson expression fitted to 
the measured mineral oil data and from the expression for constant [3 = 6.9 * 10 7 Pa, the 
same as applied in [20]. The solution from the constant [3 expression mimics the results 




Figure B.3: Pressure solutions with boundary conditions 0 W = 1.0001 and 0 e = 1. |3v = 
0.069 GPa with Sp = 12 %. r| = 0.039 Pas and u = 4.57 m/s. 
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from [20] but differs from the Dowson & Higginson expression in load carrying capacity 
with 12 %. A large difference is expected since the Dowson & Higginson expression 
represents [3 = 3.34 GPa at p p, which is significantly different from 0.069 GPa used 
in reference [20]. In Fig. B.4 the twin parabolic geometry used in [20] is shown and the 
corresponding pressure solutions are shown in Fig. B.5. The two pressure maxima are 




°0 “ * 0.0762 

x m 

Figure B.4: A twin parabolic slider with Hq = 25 pm and A = Hq m 

similar but the pressurized inlet causes the pressure maximum to the left to be slightly 
higher (1.6 kPa). There is a difference in load carrying capacity of 15 % between the two 
compressibility relations. Continuing now with the same type of geometry as in Fig. B.2 




Figure B.5: Pressure solutions across the domain with boundary conditions 0 W = 1.0001 
and 0 e = 1. [3v = 0.069 GPa with Sp = 15 %. r| = 0.039 Pas and u = 4.57 m/s. 

but with Ho = 2 pm, A = Hq and L x = 0. 1 m. 0=1 both at inlet and outlet and hence 
flooded. The corresponding dimensionless density solution where the fitted expression is 
compared with the original Dowson & Higginson expression in Fig. B.6. From left to 
right, slight liquid compression occurs and full film solutions are achieved until the gap is 
expanding to a level where the liquid does not occupy the complete gap and cavitates. The 
pressure gradient is switched off and as a consequence, the Poiseuille flow is terminated 
and the lubricant is only transported by means of Couette flow. The fractional film content 
(0) is rapidly decreasing further to the right in the domain. At the outlet boundary the 
lubricant is abruptly reformed into full film again. A significant difference in compression 
can be seen at the full film portion. Pressure solutions for the same conditions can be 
seen in Fig. B.7. The differences in pressure between the expressions are small, only 
1.7 % even though the differences in compression (the solution variable) between the to 
compressibility expressions are substantial. 

The dimensionless density solutions for the fitted and original Dowson & Higginson 
expression for a starved inlet with a fractional film content of 55 %, 0 = 0.55, can be seen 
in Fig. B.8. From the left, only Couette flow is present until the gap contracts sufficiently to 
compress the liquid and reformation into full film occurs abruptly. The reformation point 
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Figure B.6: Density solutions across the domain for the parabolic geometry with Ho=2/Ltm 
and A = Hq and boundary conditions 0 W = 1 and Q e = 1. 5g = 0.61 %, r| = 0.04 Pas and 
ii = 0.25 m/s. 




Figure B.7: Pressure solutions across the domain for the parabolic geometry with Hq =2 
y/m and A = Hq and boundary conditions 0^ = 1 and Q e = 1. Sp = 1.7 %, r| = 0.04 Pas 
and u = 0.25 m/s. 




Figure B.8: Density solutions across the domain for the parabolic geometry with Hq — 2 
and A = Hq and boundary conditions 0 W = 0.55 and 0 e = 1. 50 = 0.64 %, ri = 0.04 Pas 
and u = 0.25 m/s. 



is slightly shifted between the two compressibility expressions. Further to the right the so- 
lutions mimics the earlier case with flooded inlet conditions. The film rupture boundaries 
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appears at the same locations. In the cavitated region the equation is hyperbolic because 
of the Couette term. Hence, no information is signaled from the downstream locations 
and consequently, upstream differencing is applied. The lack of influence from the down- 
stream locations in the cavitated zones together with the the switch function give rise to 
this abrupt and somewhat unexpected reformation into full film. The corresponding pres- 
sure can be seen in Fig. B.9 where greater differences arise between the two Dowson & 
Higginson expressions, 7.4 %. Here, the expression for (3 =constant, taken for the fitted 




Figure B.9: Pressure solutions across the domain for the parabolic geometry with Hq =2 
pm and A = Hq and boundary conditions 0 W = 0.55 and 0 e = 1. ▼ : Sp = 7.4 %, ◄: 
Sp = 0.1 1 %. 11 = 0.04 Pas and u = 0.25 m/s. 

expression at 0 = 1 is added to the figure to show that the difference in solutions between 
this and the fitted Dowson & Higginson expression are negligible (overlapping lines) for 
these conditions. By inspecting Fig. B.l it is evident that the fitted expression and the 
constant (3 expression are close in a pressure range of 0-100 MPa which is reflected in 
the results where the simulations are put so that the limits of hydrodynamic pressures are 
not exceeded. In Fig. B.10 the same physical conditions are applied as in Fig. B.9 and 
solutions are obtained for different values of (3. Substantial differences in pressures are 




Figure B.l 0: Pressure solutions across the domain for the parabolic geometry with Hq=2 
pm and A = Hq and boundary conditions 0 H ’ = 0.55 and 0 e = 1 . (3 = [▼3.34/2, ◄ 3 .34 * 2] 
GPa with Sp = [8.4, 4. 8] %. r| = 0.04 Pas and u = 0.25 m/s. 

seen between the solutions indicating the importance of using a relevant value of the bulk 
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modulus for the lubricant to be simulated. 

In Fig. B. 1 1 a linearly converging film thickness is used with Hq = 1 /./m and A = 0.25 
/rm where the inlet film thickness is Hq + A and the outlet film thickness is Hq. To further 




Figure B.11: Pressure solutions across the domain for the converging gap geometry with 
Hq = 1 ^/m and A = 0.25 /j<m and boundary conditions 0 W = 0.85 and 0 C = 1 . ▼ : Sp = 22 
%, ◄: Sp = 0.54 %. Tj = 0.04 Pas and u = 0.25 m/s. 

study the effects of starvation, 0 = 0.85 is set at the inlet while ambient density is applied 
at the outlet. Also for this case the constant (3 expression produces results close to the 
mineral oil with a load difference of 0.54 %. However, the original Dowson & Higginson 
expression deviates significantly from the mineral oil with a load difference of 22%. Great 
differences are also seen in Fig. B.12 for the different values of (3. 




Figure B.12: Pressure solutions across the domain for the converging gap geometry 
with Ho = 1 A'm and A = 0.25 /jm and boundary conditions 0 W = 0.85 and Q e = 1. 
(3 = [▼3.34/2, ◄ 3.34 * 2] GPa with Sp = [25, 19] %. r| = 0.04 Pas and u = 0.25 m/s. 

The conditions at the inlet are significant in many aspects. We have seen that different 
compressibility models produce similar solutions in a pressure range reaching up to 50 
MPa but become significantly different if the inlet is starved. Let us use the fitted Dow- 
son & Higginson compressibility expression and focus on some geometrical aspects. In 
Fig. B.13 two sinusoidal waves applied on the left half of the geometry are shown. The 
inlet film thickness is Hq and Hq + A respectively but the same at the end of the two sinu- 
soidal waves. Pressure solutions for this geometry can be seen in Fig. B. 14. The geometry 
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Figure B.13: Geometries with parallel gaps, two sinusoidal waves with a difference in fre- 
quency both at (x < L/2) with Ho = 10 /jm and A = 0.5 * Hq. 



with the smaller film thickness at the inlet produces higher pressures. The results are some- 
what unexpected since the geometry with a greater film thickness at the inlet is likely to 
cause a greater compression of the fluid. In this case however, the sinusoidal geometry is 




Figure B.14: Pressure solutions across the domain for the sinusoidal geometry at (x < 
L/2) with Ho = 10 ,i/m and A = 0.4 *// 0 m and boundary conditions p w = le + 05 and 
p e = le + 05 Pa. s/ = 35 %. iq = 0.04 Pas and u = 0.25 m/s. 

able to pull a larger amount of fluid into the domain up to the point where the gap expands 
and the fluid cavitates. The non-dimensional density at the location of the first sinusoidal 
wave is shown in Fig. B.15. It can be seen that the fluid tends to get more compressed for 
geometry 2 until the earlier diverging gap leads to an earlier cavitation which affect the 
compression. 

If the same sinusoidal waves as in Fig. B.14 are squared, we get a film thickness as in 
Fig. B.16. Here the smallest film thickness is Ho- The pressure results can be viewed in 
Fig. B.17. If the inlet is flooded the geometry 1 cavitates throughout the domain producing 
zero load carrying capacity. Since the inlet film thickness for the geometry 2 is larger 
than Ho some fluid compression can occur and the pressure will rise across the domain 
producing load carrying capacity. If the geometry is pressurized, it is now possible even 
for the sin geometry to produce a pressure build-up. 
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Figure B.15: Density solutions across the domain for the sinusoidal geometry at (x < L/2) 
with flo = lOjum andA = 0.5*Ho and boundary conditions p w = le + 05 and p e = le + 05 
Pa. r| = 0.04 Pas and u = 0.25 m/s. 
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Figure B.16: Geometries with parallel gaps, two squared sinusoidal waves with a difference 
in frequency both at (x < L/2) with Hq = 10 pm and A = 0.5 * Hq. 




Figure B.1 7: Pressure solutions across the domain for the cos 2 at (x < L/2) modified step 
geometry with Ho = 10 yrm and A = 0.5 *Ho and boundary conditions p w = 2 * 10 5 and 
p e = 10 s Pa. r| = 0.04 Pas and u = 0.25 m/s. 
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B.5 Conclusions 

A generalized Reynolds equation accommodating for an arbitrary compressibility relation 
and cavitation has been derived. It has been shown that different compressibility mod- 
els can produce significantly different results in a hydrodynamic pressure range. Load 
carrying capacity is especially sensitive to the type of compressibility relation used if the 
inlet is starved. This indicates the need of employing realistic density-pressure relations. 
To produce better predictions of cavitation regions, fractional film content and pressure, 
more physical variables that are put dependent such as viscosity might be needed. Surface 
roughness close to the inlet has been shown to have important effects on load carrying 
capacity. 
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Abstract 

This paper describes a method to compute flow factors compensating for an arbitrary 
surface roughness in compressible hydrodynamic lubrication based on a homogenization 
technique. Reynolds equation is used as the governing equation and the two-scale expan- 
sion involved in the homogenization process enables the local roughness scale to be treated 
separately from the global geometry scale. With this method it is possible to compute flow 
factors for for any deterministic roughness. Measured two-dimensional surface profiles 
are used as examples. By utilizing a randomization technique it is shown that profiles hav- 
ing the same Abbot-Firestone bearing ratio also have the same flow factors, providing an 
efficient classification of surfaces in hydrodynamic two-dimensional contacts. Flow fac- 
tors are computed for the surface profiles and solutions for global problems are obtained 
and compared with the corresponding smooth solutions. 
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C.l Introduction 



It is known that the surface roughness affects the performance in various ways in a lu- 
bricated contact. Even if the surfaces are smooth on a macroscopic level, asperities on a 
microscopic level will affect the overall performance in terms of pressure build-up, film 
thickness and friction. As computer performance increases the theoretical simulations of 
lubricated contacts become more focused on treating the roughness. Since the width of a 
typical asperity is much smaller than the total length of the contact, deterministic compu- 
tations become difficult even with today’s computer power. Thus it is desirable to find a 
way to reduce the number of degrees of freedom and still treat the effects from the surface 
roughness. 

Stochastic theories such as those reported by Tzeng and Saibel [24] for one-dimensional 
transversal roughness and by Christensen [25] for longitudinal and transversal roughness 
are early approaches to deal with roughness effects in hydrodynamic lubrication. Elrod 
[26] performed a multiple-scale analysis for two-dimensional roughness where he assumed 
that the characteristic wavelength is small enough for the approach to be permitted. Patir 
and Cheng [27, 28] derived a method where the Reynolds equation is rewritten in terms 
of control volume averaged flow factors, applicable to any roughness structure. The flow 
factors are determined empirically by simulating model bearings with surfaces having cer- 
tain roughness statistics. They showed that the ratio of x and y correlation lengths plays 
an important role in characterizing the flow factors. A lot of work has been carried out 
by following the footsteps of the early flow factor approach by Patir and Cheng. More re- 
cently, variants of the Patir and Cheng flow factor computations include Lunde and To ruler 
[29] where the use of boundary conditions can be avoided, Letalleur et al. [30] where flow 
factors for sinusoidal surfaces are studied and can be retrieved analytically near contact, 
Wang et al. [31] who use the Patir and Cheng approach and incorporate a contact model 
and extend the problem into a mixed EHD model. Harp and Salant [32] who use the flow 
factor model and formulate the flow problem to treat inter-asperity cavitation. 

Recent work has paid much attention to homogenization techniques when treating the 
surface roughness. The flow problem is decoupled into two problems; the homogenized 
problem on the macroscopic scale describing the bearing geometry, and the local (or cell) 
problem on the microscopic scale describing the roughness. Examples of such work in- 
clude Kane and Bou-Said [33] who compare homogenization and direct techniques, Jai 
and Bou-Said [34] who compare homogenization and averaging techniques and Bayada et 
al. [35] who use a homogenization technique which also treats the concept of cavitation. 

In this work we present an approach to generate flow factors compensating for the 
surface roughness in the full film hydrodynamic regime using the homogenized results of 
the compressible Reynolds equation presented by Almqvist and Dasht [36]. The flow fac- 
tors can be computed for an arbitrary periodic two- or three-dimensional roughness which 
permits the use of measured surface topographies of real surfaces. In the current work two- 
dimensional surface roughness is considered and as a numerical test case, four measured 
surface profiles are used. Surface profiles are also generated artificially in a randomized 
way but with the the same Abbot-Firestone bearing ratio as the original measured surfaces. 
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The Reynolds equation is well suited for solving fluid film flows in e.g. bearings. However, 
to simulate a complete bearing including the surface roughness is not feasible because of 
the amount of degrees of freedom required to resolve the roughness. In this section a way 
to separate the effects of the surface roughness from the effects of the global geometry of 
the bearing is described. Let us assume that a surface can be quite well represented by 
using two geometry scales where one scale is much smaller than the other. With this in 
mind we can express the total film thickness as 

h(x) = ho(x) + hi(x/e), e>0 (C.l) 

where Iiq represents the global geometry scale and h\ the roughness scale which is an 
arbitrary periodic function expressed in one period with a wavelength parameter e. For this 
condition, the compressible Reynolds equation is homogenized by Almqvist and Dasht in 
[36], The homogenization process involves an asymptotic two-scale expansion for the 
solution variable and letting e — > 0. This permits the local and the global problems to 
be treated separately. Since the global problem do not involve any rapid oscillations it 
requires less numbers of degrees of freedom to be accurately resolved. 

The above discussion will be revisited later but consider now a measured two-dimensional 
surface profile r(x) with zero mean value. In the hydrodynamic contact the profile and re- 
lated measures are shown in Fig. C.l. Now the film thickness of the lubricated contact can 
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Figure C.l : Film thickness and related parameters. 



be expressed as 



h (x) = h — r. 



(C.2) 



The film thickness is constant in time and the stationary compressible Reynolds equation 
constitutes the basis in this work, viz. 



/ p/z 3 dp\ u d (p h) 

\12r|r/A'/ 2 dx 



d_ 

dx 



(C.3) 
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where p is pressure, T| viscosity and p is the density. It is convenient to express the equation 
in terms of density as dependent variable so we define the dimensionless density as 

0 = — (C.4) 

Pc 

where p r is the density at ambient pressure p = p c . We assume that p = p(Q) and hence it 
is possible to apply the chain rule on the pressure gradient, i.e.. 



dp dp dQ 

dx dQ dx 



(C.5) 



Substituting the expression for the pressure gradient and 9 into the Reynolds equation 
yields: 

d f Qh? dp dQ\ ud(Qh) 
dx y 12r| dQ dx ) 2 dx 

Assuming further that the bulk modulus, defined as 



°dQ’ 



(C.7) 



can be treated as a constant yields the following expressions relating pressure to density: 



P = Pc + pln(0) 



(C.8) 



and the equivalent expression relating density to pressure: 



0 = exp 



P - Pc 

P 



(C.9) 



Note that the constant bulk modulus gives rise to an exponential expression for Q(p) which 
is an acceptable compressibility relation for a limited pressure range, see [37] for a more 
comprehensive discussion regarding this matter. Under this condition Eq. (C.6) can be 
written: 



d / 3 dQ\ _ A d{Qh) 
dx \ dx J dx 



(C.10) 



where A = 6 t|m/(3. To obtain a homogenized equation in a more general form independent 
of A we define the following dimensionless parameters: 



H = 



h 
Ji ’ 




r max 



r max _ 6r| UX _ p 

—=— ; x = — - 77 ; p = — . 

h P ymax) Pc 



(C.ll) 



Also let r = r(^) which now corresponds to h\ in Eq. (C.l) where q = x/e. Dividing 
Eq. (C.2) by h and using the above expressions the dimensionless film thickness equation 
can be written as 

ff a (5) = l+a(G(5)-l), 0 < a < 1, (C.12) 

where the subscript indicates that a is a parameter. As a goes from zero to one, the 
clearance h goes from zero to infinity. It should be noted that the mean value of H always 
equals unity and as the a decreases, H becomes smoother and equals to unity for a = 0. 
Substituting H a , p and x into Eq. (C.10) the dimensionless form of the equation becomes: 



d_ 

dx 




d{QH a ) 



dx 



(C.13) 
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The corresponding homogenized equation (see Almqvist and Dasht [36]) yields: 



d 

dx 





(C.14) 



where 

<!>*« = J H a (5) (l + d\ and 

(C.15) 

\)/ and x are solutions to the local problems, i.e. the local roughness scales which are given 



by 



d_ 




3 (^) 



and 



— ( ft 3 _ dH a 

dt; V “ ^ J - dt ; ' 



(C.16) 



The dimensionless Reynolds equation has lead to the local problems being independent 
of A. Since the local problems are solved at a unit cell domain, physical spatial properties 
of the roughness are diminished which explains the choice of non-dimensionalization of x 
in Eq. (C. 11). Simply speaking, the influence of x disappears in the solutions of the local 
problems. 

The dimensionless film thickness Eq. (C.12) has led to the homogenized equation be- 
ing independent of x but dependent on the parameter a. This permits the local problems to 
be solved for 0 < a < 1 producing flow factors (|) A and § s as functions of a which add the 
contribution from different roughness G(^) to a global coarse homogenized problem. With 
a database of previously computed flow factor the non-dimensional homogenized equation 
can be directly solved on the non-dimensional coordinate x to achieve the solution 8. In- 
stead of using the somewhat inconvenient definition of x the dimensionless homogenized 
equation can be directly scaled into dimensional form 



= f (♦.W*WA0), 



(C.17) 



where (j)* and § s are retrieved from the parametrized § X(J and <\> Sa by e.g. table look-up to 
become functions of x. This equation could be compared with the one derived by Patir and 
Cheng [28] but with the use of different flow factors § x and (]) s . 



C.3 Results 

C.3.1 The Effect of £ 

In the two-scale approach involved in the homogenization process the wavelength measure 
8 of the local roughness scale is set to approach zero. This would correspond to a deter- 
ministic solution for a geometry with a periodic roughness wavelength approaching zero. 
The thin film approximation involved when solving the Reynolds equation only permits 



60 



PAPER C. 



solving a problem with a large width to height roughness ratio. However, to verify the ho- 
mogenized solution, a corresponding deterministic solution could be obtained where the 
roughness wavelength decreases to a finite level. The deterministic solution should then 
approach the homogenized. Consider the film thickness equation 

h = ho +Asin (C.18) 

where Iiq is a convergent gap and A = min(/zo) / 10. This film thickness is shown in Fig. C.2 
with sinusoidal waves of different wavelengths. Deterministic and homogenized pressure 
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Figure C.2: A convergent gap geometry model with sinusoidal roughness with e = 2 5 
and 2“ 3 . 



solutions for this geometry are shown in Fig. C.3. It can be seen that the deterministic 




(a) Full domain (b) Magnified area 

Figure C.3: Pressure distribution along a converging gap with sinusoidal roughness for the 
homogenized and deterministic solutions. 

solutions approach the homogenized as the wavelength decreases. Good correspondence 
between the deterministic and the homogenized solutions is achieved for a finite value of 
the wavelength. This motivates the treatment of the limit e — ► 0 in the homogenization 
process. 



C.3.2 Flow Factors 

In this section four measured two-dimensional surface topographies are used in the numer- 
ical examples. The machining process applied to the surfaces and some of their surface 
parameters are given in Table C.l. It can be seen that surface A differs a lot from the 
other surfaces in R a value. Surfaces B and D are quite similar but differ significantly in 
skewness. The Abbot-Firestone bearing ratio for the four profiles can be seen to the left in 
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Table C.1 : Data for the four original surface profiles. 



# 


Machining 


Ra 


R, 


SK 


K 


A 


Unidirectional ground 
+ Phosphated and dephosphated 


0.93 


5.0 


-0.99 


4.3 


B 


Unidirectional ground 


0.29 


1.2 


-0.49 


3.3 


C 


Unidirectional ground 
+ 3 x chemically deburred 


0.13 


0.58 


-0.010 


3.6 


D 


Unidirectional ground 
+ shot-peened 


0.37 


1.6 


0.032 


3.0 





Figure C.4: Abbot-Firestone bearing ratio for the surfaces in dimension (left) and dimen- 
sionless (right). 



Fig. C.4. This cumulative height distribution contains all height information about a sur- 
face such as for example R a , R q , skewness and kurtosis but no spatial information. Since 
the simulations will consider the dimensionless surface roughness we show the cumula- 
tive height distribution in terms of G, for clarity, to the right in Fig. C.4 where the mean 
value of G equals unity. With this dimensionless representation the absolute height values 
are insignificant and only the shape of the Abbot-Firestone bearing ratio is relevant. The 
results to be discussed are based on these dimensionless height distributions. 

The surface measurements consist of 512 points equally spaced and distributed on 
a length of 4.2 x 10~ 4 m. To compute the flow factors, see Eq. (C.15), the number of 
degrees of freedom of the measurement has to be increased. The computational grid has 
to be sufficiently dense in order to resolve the roughness on a cell. 

For this reason extra grid nodes are linearly interpolated between the points of the sur- 
face measurements. Fig. C.5. Because of the linear interpolation no information is added 
or removed and the surface representation will become identical to that of the the mea- 
surements, only more degrees of freedom are added to the solution process. The demand 
for more degrees of freedom is illustrated in Fig. C.6(a) where the flow factors § s are 
shown for surface A as functions of the parameter a for increasing number of grid nodes 
N. Solutions from 50 values of a, equally spaced between a = 10~ 9 and a = 1 — 10~ 9 are 
used 1 . The dashed line represents the original profile and it is seen that the flow factors 

1 Mechanical contact occurs for a = 1 
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Figure C.5: A small part of the original surface and linear interpolation with one point. 





(a) Surface A (b) Surface B 

Figure C.6: Flow factor with increasing number of interpolation points where dashed line 
represents the original surface with 2 9 points and highest resolution \sN = 2 16 . A total of 8 
curves are shown in each plot 



become more and more sensitive to grid resolution as a approaches unity. Because of the 
way that H is non-dimensionalized, the mean value of H is always unity and the shape of 
H becomes smoother as a approach zero which means that the demand on resolution de- 
creases. Particularly, at the extreme of a = 0 the local problem can be solved analytically. 
The same conditions are applied for surface B in Fig. C.6(b). It is clear that this profile is 
less sensitive to grid resolution than surfaced. To display the change of the flow factors as 
the grid spacing is reduced, we define the difference in solutions between two consecutive 
grid spacings 8 and and its double spacing 28 respectively as 




(C.19) 



where (f» is taken at some value of a. A visualization for convergence measure is seen in 
Fig. C.7. The linear behavior for the convergence in flow factor for a = 0.49 at higher 
resolutions indicates that the error reduction is in agreement with the effective order of the 
discretization scheme which is not that obvious for a close to unity. However, the flow 
factors for a film thickness approaching zero are somewhat uncertain regardless resolution 
since the hydrodynamic regime is to some extent violated. It can also be seen that the 
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(a) a ~ 0.49 



(b) a ~ 0.49 



10 ° 



<-£ 1(T 5 



ter 10 

2 




(c) a ~ 1 (d) a ~ 1 

Figure C.7: Flow factor convergence as a function of spatial resolution. 



convergence rate for surface A is generally smaller than for the other surfaces. 



C.3.3 Surface Classification 

It is desirable to find classes of surface roughness that affect the pressure in a hydrody- 
namic contact similarly, i.e. having the same flow factors. In the Abbot-Firestone bearing 
ratio all surface height parameters are known but no information about spatial parameters 
is revealed. In this work the Abbot-Firestone bearing ratio is Considered as a suitable tool 
for surface classification. To do this the Abbot-Firestone bearing ratios from the origi- 
nal surface profiles are used and the points at discrete levels of the curves are spatially 
distributed in a random way to produce surfaces having the same Abbot-Firestone bear- 
ing ratio. In order not to violate the thin film approximation, the distribution of points is 
controlled in certain ways. A fully randomized surface could have the have highest and 
lowest points as neighbors. This is accomplished by introducing some rules so that points 
are allowed to be positioned spatially close to existent points at a higher height level. If 
that is not possible, new peaks are formed from where points can be positioned at a lower 
height level. Hence, the points are not fully randomly distributed. This is because the pro- 
files would violate the thin film approximation, with the possibility of highest and lowest 
points being neighbors. Flow factors for the original and five random surfaces are shown 
in Fig. C.8. At the top two figures the original measured profiles with N = 2 9 are used. The 
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Figure C.8: Flow factors against a for all surfaces. The resolution \s N = 2 9 top and 
N = 2 16 bottom. For each original surface (dashed), five randomizations of that surface are 
also plotted. 



randomized profiles are fairly well clustered together but the original profiles differ signif- 
icantly. As discussed before, more degrees of freedom than those of the original profiles 
need to be introduced, especially for a close to unity. At the bottom two figures interpo- 
lation is applied where N = 2 16 . All surfaces having the same Abbot-Firestone bearing 
ratio coincide and it is clear that for a sufficiently resolved profile, the Abbot-Firestone 
bearing ratio is by itself determining the result. Hence the Abbot-Firestone bearing ratio 
holds the information necessary to predict the influence of surface roughness on the flow 
factors introduced by the homogenization process. 

To get a better idea of the deviation between the 6 profiles, i.e., the five randomized 
profiles and the original profile, the standard deviation of § x is computed for each value 
of a. Fig. C.9. The standard deviation for each surface reaches a maximum at a close 
to one. The same applies for standard deviations of (}>< but with a maximum of less than 
1.2 x 1 0 6 . The maximum standard deviation decreases as N increases which is obvious 
from Fig. C.10. This confirms that the Abbot-Firestone bearing ratio acts as a good clas- 
sification. To confirm that the maximum standard deviation is small enough, the homog- 
enized equation is solved for a dimensional converging gap geometry where the smallest 
film thickness is represented by a = 0.99 and the largest film thickness by a = 0.6 for 
all four surfaces. This is in order to cover a range where the standard deviations of flow 
factors are the largest for all four surfaces, see Fig. C.9. The load carrying capacity is then 
used to monitor the total error. Each surface ( A,B,C,D ) is interpolated to N = 2 16 grid 
nodes and the standard deviation for load carrying capacity was computed. For all sur- 
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a 



Figure C.9: Standard deviations of (|>. T for the random and original profiles where N = 2 16 . 




N 

Figure C.10: Maximum standard deviations of ^ in the a-range. 



faces, the value of this measure was found to be less than 0. 1 % of the mean load carrying 
capacity. 



C.3.4 Application Examples 

Some dimensional bearing examples will be shown where the computed flow factors are 
used to solve the homogenized problem. The non-dimensional flow factors for the original 
profiles are used and scaled into dimensional ones (see Section C.2) allowing the homog- 
enized equation to be solved for the dimensional spatial domain. The flow factors needed 
for the corresponding film thickness are discretely retrieved by cubic spline interpolation 
from the computed flow factors. The film thickness notations are in terms of the mean 
nominal separation h. In all simulations, the domain length in the jr-direction L x = 0.1 m, 
p = 0.14 Pas, (3 = 10 9 GPa, u = 1 m/s, p c = 10 5 Pa which also are the inlet and outlet 
Dirichlet boundary conditions. It should be noted that side leakage is not considered since 
the one-dimensional Reynolds equation is used and therefor, the load carrying capacity 
becomes greater than for two-dimensional problems. 

To verify if the local problems are sufficiently resolved in terms of N to produce suf- 
ficiently accurate solutions to the dimensional homogenized problems we compute the 
convergence in the load carrying capacity between two consecutive grid spacings 8 and its 
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double spacing 28 respectively as 



W 6 



W 5 _ W 25 

P 



where the load carrying capacity is computed as 



W = f p dx. 
Jo 



(C.20) 



(C.21) 



The load convergence for the converging cap with min(/z) = 5 pm and ma x(h) = 6 pm is 
shown in Fig. C.ll from where it can be seen that the critical limit is IV = 2 11 and the 




Figure C.1 1 : Load convergence as a function of grid resolution for a converging gap. W < 

1 .8 x 10~ 5 for all surfaces at highest resolution with N = 2 16 . 

change at the highest resolution is considered as sufficiently small. Thus, N = 2 16 is the 
resolution used when displaying all dimensional results. 

The minimum global film thickness min (/io(x)) is used as a solution variable, deter- 
mined in order to achieve force balance for a given load. 

Consider a converging gap geometry with max (ho(x)) = min(/ro(x)) + 1 pm. In Fig. C.12 
the minimum value of Iiq is plotted as a function of load carrying capacity for the rough 
profiles and a smooth, denoted S. The ordering of the minimum film thickness follows the 
ordering of R a values for the roughness profiles, with surface A standing for the smallest 
film thicknesses and the smooth S for the largest, respectively. It is also clear from the 
figure that as the film thickness increases the influence of roughness decreases. 

Now consider a parabolic geometry shown in Fig. C. 13. Since atmospheric pressure is 
applied at the boundaries, cavitation is likely to occur at some point in the diverging part. 
For this case we use the approach of treating cavitation by Vijayaraghavan and Keith [20]. 
They use a Reynolds equation in a different form than presented here. This means that 
the homogenization process does not apply for that equation. However, since the equation 
only differs in a switch function applied to the Couette term the equation is identical to 
the one which is homogenized in the flooded region where the switch function equals one. 

If Eq. (C.17) is modified by the switch function to treat cavitation, the global problem 
will treat cavitation but not the local problem. This leads to vanishing of § x at cavitation 
but the influence from <j> 5 will remain, pressure distribution obtained from the modified 
homogenized equation for surface A is given in Fig. C.14. Pressure is rising when the 
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Figure C.12: The minimal film thickness as a function of load carrying capacity for the 
converging gap geometry. 




o “ * o. 1 



x m 

Figure C.13: A parabolic geometry with max(/io) =min(/ro) + l x 100 ^m. 




Figure C.14: Pressure across the domain for the parabolic geometry when surface A is 
used. W = 100 kN/m. 



gap is converging and in the diverging part where x « 0.07 m cavitation occurs. The 
effect of surface roughness on film thickness Iiq as a function of load carrying capacity is 
shown in Fig. C.15. As for the converging gap geometry the rough surfaces give smaller 
film thickness Iiq compared to the smooth geometry. And the same order in which rough 
surfaces affect Iiq applies here too. 



C.4 Conclusions 

In this paper we consider a homogenization technique for the compressible Reynolds equa- 
tion. The approach involves an asymptotic two-scale expansion of the film thickness. For 



68 



PAPER C. 




Figure C.15: The minimal film thickness as a function of load carrying capacity for the 
parabolic geometry. 



a lubrication problem involving rough surfaces, this method provides a mathematical ap- 
proach that separates the hydrodynamic effects at the local scale, representing the surface 
roughness, from the hydrodynamic effects at the global scale, representing the bearing 
geometry. The conclusions are as follows: 

• A novel method of computing dimensionless flow factors for arbitrary surface rough- 
ness using a homogenization technique has been developed. 

• The method provides flow factors for arbitrary global geometries within the full film 
regime. 

• It has been shown that different profiles having the same Abbot-Firestone bearing 
ratio give the same flow factors, thus providing an efficient classification system for 
surfaces intended to operate in full film lubrication. 
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